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Abstract 

We investigate variational methods for finding approximate solutions to the 
Fokker-Planck equation, especially in cases lacking detailed balance. These 
schemes fall into two classes: those in which a Hermitian operator is con- 
structed from the (non-Hermitian) Fokker-Planck operator, and those which 
are based on soluble approximations to this operator. The various approaches 
are first tested on a simple quantum-mechanical problem and then applied to 
a toy Fokker-Planck equation. The problem of a particle moving in a po- 
tential and subject to external non- white noise is then investigated using the 

formalism developed earlier on in the paper. 
PACS numbers: 05.40.+j, 02.50.Ey, 03.65.Db 
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I. INTRODUCTION 



Variational schemes, such as the Rayleigh-Ritz procedure Jlj] are well known methods 
for finding approximate solutions to the Schrodinger equation. There has been far less 
work done in applying analogous procedures to the Fokker-Planck equation, even though 
the equations share a similar structure. One of the reasons for this situation is that the 
differential operator in the Fokker-Planck equation is not, in general, self-adjoint, making 
the formalism more complicated: eigenvalues are not necessarily real, nor are the right and 
left eigenfunctions equal. Another reason is the existence of a zero eigenvalue, corresponding 
to the steady-state of the system. Under certain conditions (the "potential conditions" ) the 
determination of the steady-state probability distribution function reduces to quadratures. 
Since all eigenvalues have a non-negative real part 0, in this case the lowest eigenvalue and 
eigenfunction are known exactly, and approximation schemes have only to be developed for 
"excited" states. The formalism for this has been developed and applied to a number of 
problems |||| . However, for many systems of interest the potential conditions do not hold, 
and as a consequence the steady-state distribution cannot be determined in closed form. It 
is desirable to have variational schemes in this case as well. Some work has been carried out 
by Seybold [|J, but no systematic discussion exists in the literature. The purpose of this 
paper is to investigate the usefulness of a number of variational schemes in the case where 
detailed balance does not hold. These are illustrated on simple two-dimensional systems, 
which are the simplest systems for which potential conditions need not hold. 

The outline of the paper is as follows. In section two we review some general formalism 
for the Fokker-Planck equation and in addition discuss the potential conditions and show 
how they lead to an exact solution for the stationary probability distribution. Section 
three introduces the variational methods that are the heart of this work In section four, 
we consider some examples of these schemes in action, including their application to the 
colored-noise problem. 
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II. FORMALISM 



In this section we first present some of the formalism associated with the Fokker-Planck 
equation and then treat the potential conditions, demonstrating how they lead to an exact 
solution for the steady-state probability distribution. Here, we will see that making the 
Fokker-Planck equation look as much as possible like quantum mechanics is very natural. 

A. General Formalism 

The Fokker-Planck equation for a system with N degrees of freedom is 

d t p(x,t) = Lp(x,t), (2.1) 

where L is a differential operator of the form: 

L = -dMx) + didjBijix), (2.2) 

where the summation convention is understood. Here Ai(x) and B^{x) = Bji(x) are the 
drift vector and the diffusion matrix respectively, d\ means d/dxi and i,j = 1,2,...,N. 
Looking for separable solutions of the form: 

p{x,t) = P n (x) exp{-A n t} (2.3) 

leads to the eigenfunction equation: 

LP n {x) = -\ n P n (x), (2.4) 

where we have assumed for convenience a discrete eigenvalue spectrum with n = 0, 1, 2, .... 
Since L cannot, in general, be brought into Hermitian form, the eigenfunctions Q m (x) of 
the adjoint operator have to be found: 

L ] Q m (x) = - \ m Q m (x), (2.5) 
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as well. The set of eigenvalues in ( |2.4|) and ( |2.5| ) are equal; moreover, since the operator 



is real, the complex conjugate of an eigenvalue is also an eigenvalue. The set of functions 
P n (x) and Q m (x) are bi-orthogonal: 

J d N x P n {x) Q m {x) = 5 n>m . (2.6) 

Equation ( |2.3| ) shows that the stationary state corresponds to Ao = i.e. Po(x) = Pst(x)- 

From the form of L we see that Qo(x) is a constant, and then from ( |2.6| ) that Qo(x) = 1. 
Henceforth, we will frequently denote operations such as those found in eqs. Q2.4|), ( |2.5| ) and 
( p.6| ) in the following notation: 

L \P n ) — \Pn) 
(Qm\ L X m (Q m \ 

(Qm\Pn) 3n,m- (2.T) 



B. The Potential Conditions 

The study of the Fokker-Planck equation simplifies considerably if the drift vector Ai(x) 
and the diffusion matrix P>ij(x) satisfy the so-called potential conditions mentioned in Section 
I. There are various ways to introduce these conditions, but perhaps the simplest is first to 
note that the Fokker-Planck equation may be written in the form of a continuity equation: 
dtp(x,t) + diJi(x,t) = 0, where Ji is: 

Ji{x) = [Ai(x) -djBijix^p&t) (2.8) 

and called the probability current. Stationarity implies V ■ J_ = 0, but if it is also true that 
J = 0, then 

dj [Bij(x) p s t(x)] = Ai(x) p s t(x). (2.9) 

(It turns out that an equivalent statement is P n (x) = p s t{x)Q n (x) for all n.) If in addition 
the diffusion matrix B has an inverse then 
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d(x) = di\np st (x) = B^\x)[Aj(x)-d k B jk (x)] (2.10) 

is a gradient. Necessary and sufficient conditions for this are: 

^ = 22, (2.11) 

with i,j = 1, ...,7V. Provided these conditions hold, the determination of the stationary 
probability distribution is reduced to quadratures: 

Pst(x) = expjy~d^ Ci(x!_)\. (2.12) 
Furthermore when fl2.9p holds, the operator 

C = ( Ps t(x)y 1/2 L (p st (x)) 1/2 (2.13) 
is Hermitian M with the real eigenfunctions 

Vfe(30 = : P "^ 1/2 = (Pst(x)) 1/2 Q n (x), (2.14) 

and the eigenvalues are non-negative. 

Of course, these simplifications occur only in exceptional cases; for most cases the oper- 
ator L cannot be brought into Hermitian form as in ( |2.13j ). Nevertheless, it is still useful 
to define the operator C, since it turns out that splitting this operator (as opposed to L) 
up into Hermitian and ant i- Hermitian parts is natural in the sense that it is dictated by 
the underlying temporal symmetries of the system 0. This decomposition can be achieved 
explicitly by writing Ai = Af + Aj~, where Af is defined by Af = p^djlB^p^]. Denoting 
the operator of the form fl2.2|) , with A replaced by A + , L by L + and defining L~ = —diA~, 
the original Fokker-Planck operator may be written as L = L + + L~ with the corresponding 
£ operators defined by ( |2.13j ) having a similar decomposition: C = C + + £~ . By construc- 
tion C + is Hermitian, since A + satisfies the conditions ( |2.9[ ). It is also easy to check that 
C~ is ant i- Hermitian using L~p st (x) = 0. 

In general, the operators £, and will have a different set of (complex) eigenfunctions: 
from (2.4), (|2.5| ) and ( 2.13|) the eigenfunctions of these operators are: 




and Xn(x) 



(Pst(x)) 1/2 Q n (x). 



(2.15) 



(Pstfe)) 



1/2 



From (|2.6| ) they satisfy the ort honor mality condition 



(Xm\tpn) = / d X1p n (x) Xm(l) = 6, 



(2.16) 



The eigenvalues, although complex, have a non-negative real part as expected on phys- 
ical grounds. The eigenfunction of both C and C) corresponding to the zero eigenvalue is 



Variational approaches are quite useful for finding ground- state energies in quantum me- 
chanics because: i) an error of order e in the variational wave function results in an error of 
order e 2 in the variational energy; and ii) the true ground-state energy is known to be lower 
than the variational energy. Since the lowest eigenvalue of a Fokker-Planck operator is iden- 
tically zero, the focus becomes the eigenstates (the lowest being the stationary probability 
distribution) and/or excited eigenvalues (the first excited eigenvalue being the reciprocal of 
the longest relaxation time of the system H). 

Some work on Fokker-Planck variational approaches has been carried out, mainly on 
one-component problems. In such cases the potential conditions automatically hold, the 
stationary probability distribution is known exactly, and the interest was naturally centered 
on higher eigenvalues. In these circumstances, progress is made by obtaining the Hermitian 
operator £ and proceeding as in quantum mechanics. In calculations of the lowest non-zero 
eigenvalue, the nice features of a variational approach can be recovered simply by ensuring 
that the variational state is orthogonal to the stationary probability distribution. However, 
since most problems do not obey the potential conditions, variational schemes applicable to 
these more general situations need to be developed. 



(P*t(x)) 



III. VARIATIONAL APPROACHES 
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A. Constructing a Hermitian Operator: The L'L Approach 



When the potential conditions do not hold, one can still mimic quantum mechanics by 
constructing a self-adjoint operator from the Fokker-Planck operator L. A natural choice 
is since it shares right eigenstates with L and its eigenvalues A n are related to those 
of L through A n = A* A n . With a Hermitian operator in hand one simply proceeds as in 
quantum mechanics. In addition to the drawback that there is no longer any interest in 
finding the lowest eigenvalue (it is identically zero), another disadvantage of this scheme 
compared to Rayleigh-Ritz in quantum mechanics, is that the calculation requires twice 
as many operations — first operating with L and then with L\ rather than just with the 
Hamiltonian H. 

It is natural to consider minimizing (Pq t \L^ ' L\P^) / '(Pq t \Pq t ) Q, where Pq t (x) is a trial 
distribution, as a means for finding a reasonable approximation to the steady-state distri- 
bution. This quantity is bounded from below by zero and equals zero only if the choice 
Pq(x) is equal to the exact steady-state distribution Po(x). 

Note that when examining the eigenvalue equation LP n = —X n P n one has the freedom to 
transform the equation to LP n = — A n P„ where L = R~ X I 2 LR}I 2 and P n = R~ 1 ^ 2 P n provided 
R has no zeros. Hence, one alternative would be to minimize (Po r |L^L|PQ r )/(PQ r |PQ r ). 
Calculational convenience sometimes dictates using R = P^. 

The first excited eigenvalue Ai (the reciprocal of the longest relaxation time of the system) 
is accessible to a straightforward variational approach provided |Pi) has a different symmetry 
than |P ). In such cases (P* r |L'* L\P^) / (P 1 tr |P 1 tr ) yields an upper bound on Ai provided one 
can guarantee that |Pi r ) is orthogonal to \Pq). 

It is possible to lower the bounds on A n by improving |P^ r ) (making a better choice or 
one with more variational parameters). A systematic procedure, which is similar to Langzos 
tridiagonalization 0, uses: 

\V tr ) = |P tr ) + aL\P tT ) + (3L 2 \P tr ) + . . . , (3.1) 

as the improved trial distribution, where a, j3, ... are additional variational parameters. Fi- 
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nally we note that it is possible to get some sense of how good the variational estimate is, 
since in the Fokker-Planck problem the lowest eigenvalue should be zero. In addition, for 
the higher eigenvalues there exists a procedure which finds lower as well as upper bounds 
on the eigenvalues |3]||. 



The other variational approaches we will discuss have their foundations in perturbation 
theory. If the Fokker-Planck operator at hand is in some sense close to one which is 
solvable, then a perturbation approach is a viable scheme of approximation. One can make 
a perturbative expansion variational by introducing parameters into the solvable part and 
eventually choosing them to minimize some measure of the "nearness" of the two problems. 

The perturbative approach begins by splitting the Fokker-Planck operator as follows: 



where is the operator of some completely solvable problem, i.e., all of its eigenvalues A^ -* 
as well as its right and left eigenf unctions, P^(x) and Q$(x) are known. The variational 
aspect of this approach requires 

L (o) 

to contain some as yet undetermined parameters. Note 
that e need not be small and is used here primarily as a counting device. Expanding the 
eigenfunctions and eigenvalues of L as a series in e yields: 



(l<°> + eL«) (|Pf> + C |P«> + ...) = - (A?> + + . . .) (|Pf ) + e\P^) + ...). 



B. Variation within Perturbation 



L = L< > + eL« 



(3.2) 



(3.3) 



Manipulating the terms of order e is the usual ways leads to: 




e (Q^\Pi 0) ) 

v \p^)(Q^\lW\p^) 

L> (\(0) x(0U 




(3.4) 



(3.5) 



(Q ( n 



(3.6) 



s 



The terms of order e 2 lead to: 



K - e 2. 77(5) (o), > ( 3 - 7 ) 



and so on. 

Notice that the structure of the perturbation expansion is such that Ao = and (Q \ = 1 
at every order, as these results are exact. Furthermore, it is worth remarking that the 
expansion for \P ) does not necessarily remain everywhere positive order- by-order; thus 
there can arise difficulties in interpreting a truncated expansion of \P ) as a probability 
distribution. Another feature of the perturbation expansion to note is that if the A^'s are 
all real, then the perturbation expansion for the eigenvalues remains real (provided there are 
no degenerate A^'s which would necessitate degenerate perturbation theory); that is, the 
imaginary part of the eigenvalue is inaccessible to a perturbation theory that begins with 
purely real eigenvalues. 

If L has been chosen so as to satisfy the potential conditions, then it is convenient to 
consider the perturbative expansion of the operator transformed by Pq°^: 

(^ + dU)(m + e\^) + ...) = -(Ai ) + e A« + ...)(|^))+ e |^) + ...) (3.8) 

where 

e<»> = (p„ < ° , (£)) _1/2 l ,0, (p„ <0, (£)) 1/2 

= " l - J 1/2 . (3.9) 

(Po (0) (-)) / 

This is because the left and right eigenstates at zeroth order are identical, which is of 
considerable calculational convenience. 

Of course, the more natural rotation is by the true stationary probability distribution, but 
that is unknown. In the perturbative approach it is possible to perform this transformation 
order-by-order as follows: 
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£ = (n (0) + eP (1) + • • •) ~ 1/2 {L {0) + eL«) (P (0) + eP (1) + . . .) 



1/2 



(3.10) 



where 



£ (o) = (p Q (o))- 1/2 L (o)(p Q (o)) 



1/2 



P { 



(o)\- 1/2 



L (o) 



(1) 



-in 



(0)\ 1/2 



(3.11) 



(3.12) 



and so on. Then one could proceed with 
(C®+dCF>+..) (|< 0) ) + e|*W) + ... 



+ eAW + ...)(|*f)+ e |^) + 



(3.13) 



This approach is rather cumbersome, especially when it is noted that finding requires 
knowing P^ . 

Since the lowest eigenvalue in the Fokker-Planck perturbation expansion is zero order-by- 
order, some other quantity must be found to vary. It is desirable, though not necessary, to 
vary a bounded expression, as the bound helps to ensure a sensible result. This motivation 
led us to consider varying (Pq^IPq 1 ^) or alternatively (\&o l^o ) which is positive semi- 
definite by construction. We call this approach the "minimal corrected wave function" 
criterion or MCW. The variation seeks to minimize (P (1) |P (1) ) - - the philosophy beine; 
that the closer \Pq°^} is to the actual p s t(x), the smaller its correction should be. The 
procedure might be continued by minimizing (P |P ) to determine the parameters to be 
inserted in \Pq°^) + |Po ), and so on. To calculate Ai, (P^lPi""^) could be minimized to fix 
the parameters to be used in A^ + Ai + ... + \[ . 

Another method based on perturbation theory with undetermined parameters has been 
used by Edwards and co-workers on problems such as polymers with excluded volume || 
and the Fokker-Planck formulation of the KPZ equation fllQfl . Stevenson has dubbed it 



the "fastest apparent convergence" criterion or FAC [|TTJ. First a number of terms in a 
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perturbative expansion of the quantity of interest are calculated. FAC then assumes that 
the zeroth-order term (which depends on the input parameters) is exact; and so the rest 
of the expansion is set to zero. This last step determines the unknown parameters to be 
substituted into the zeroth-order term. Note that FAC is not a variational approach as the 
parameters are not determined by varying. 

A scheme with the same starting point which is variational is the so-called "principle 
of minimal sensitivity" or PMS [11]. After obtaining a truncated perturbation expansion, 



one varies it with respect to the undetermined parameter (s). Note that these parameters 
were introduced artificially and that the actual answer should not depend on them; however, 
any truncated expansion does depend on them. The PMS philosophy is then to search for 
the result that is "least sensitive" to the parameters - and hence the variation. In a few 
select simple examples the PMS procedure has been proven to yield a convergent series of 



approximations even when the underlying perturbation expansion is asymptotic ||12|| , but 
the general conditions for which it does so remains an open problem. 

C. Comparison of Approaches on An Example from Quantum Mechanics 

Let us test these approaches on a well-known problem from quantum mechanics, the quar- 
tic oscillator: H = — ^d 2 /dx 2 + |x 4 using the harmonic oscillator H^°> = — ^d 2 /dx 2 + ^u 2 x 2 
as the basis for the perturbation expansion. By dimensional arguments, the eigenvalues of 
H can be seen to be proportional to g 1 ^ 3 , and hereafter we scale this factor out. Let us focus 
our attention on the ground-state energy; the result is known to be E$ irect = 0.420805 .... 
Applying the MWC, FAC and PMS [|I~2"] , |l"3|l approaches to standard Rayleigh-Schrodinger 



perturbation theory (calculated to third order for the energy and second order for wave 
function) yield the results seen in Table 1. We also include the Rayleigh-Ritz result and its 
first Langzos-like correction (RRL). 

Let us make a few observations. PMS and RRL are identical at first order. At second- 
order the FAC and PMS approaches have no physical solutions; in PMS one can search for 
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inflection points when no extrema are found — at second order this yields 0.42143. Notice 
that not only are the variational approaches better than FAC at first order but also they 
were improved by going to second or third order, while FAC got worse. To the orders 
calculated here, PMS has led to the best results. 



IV. SOME FOKKER-PLANCK EXAMPLES 

In this section we will consider a few examples of the variational schemes applied to 
Fokker-Planck problems. First, we will apply the techniques to a toy model. Then we will 
examine a more difficult problem seen recently in the physics literature — the colored-noise 
problem. 



A. A Toy Model 

The toy model we consider is the Fokker-Planck equation corresponding to the following 
two coupled, non-linear Langevin equations: 

— = -uxi - gx x x 2 + 771(f) 

^ = -vx 2 + V2 (t), (4.1) 
where the noise is Gaussian-correlated with zero mean and the following correlation: 

(i*(f)ifc(0> = 2D5 l3 5{t-t>). (4.2) 
The associated Fokker-Planck operator is: 

/ <9 2 d 2 \ ( d d \ J d\ 

*** = D fe + ^j + H 2+ ^ + ^j + ^l 1+ ^j- (4 - 3) 

The first thing to note is that it does not satisfy the potential conditions. Since they 
are not satisfied, let us apply each of the variational methods suggested above to obtain 
the steady-state values of (x\) and For the approaches based on perturbation theory 

we will need an operator for which the eigenvalue problem is solvable. We choose as L^oy 
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the Fokker-Planck operator for an Ornstein-Uhlenbeck process (the harmonic oscillator of 
Fokker-Planck problems): 



d d I d d 

L{% = vd l — I + ud 2 — I + v[2 + x l — + x 2 —). (4.4) 
ox{ OX2 \ OXi OX2 



The eigenstates of L^ y are: 

""" aW ~ (2».+«+%=n 1 ! n . 2 !<i 1 <i 2 ) 1 / 2 P l 2d, 2rf 2 /' i4 '" J 
and its eigenvalues are A n i >n 2 = — (ni + 7^2)^ where ni, n 2 = 0, 1, 2, .... As we are interested 
in spacial quantities, we have selected parameters di which affect the spatial distribution 
Po°\x) but leave the eigenvalues of L[^ y unchanged. 

The details of the variational calculations we performed on the toy model can be found 
in the appendix. The results are shown in Table 2 along with the results from a simulation 
of the toy-model Langevin equation (|4.1| ). 

The simulation algorithm employed a second-order Runge-Kutta method to evolve the 



equations. [0 After a sufficiently long evolution, dependence on the initial conditions is 
lost. We simulated the equation for 1 000 000 realizations and extracted the averages 
(xf(tf)) where Xi(tf) is the final position of the simulation. The results of a simulation with 
D = 0.5, v = 1.0 and g = 3.0 are (x\) = 0.2652280 and = 0.5003104. 

Note that the values for d\ and d 2 found from varying ('^o^lV'o 1 '') were substituted into eq. 
( |A9|) at first order. The methods based on perturbation theory all agree that (x^) = 0.5 and 
are better than the (^o°' ) |4oy4oy|V ; o ' ) ) variation on this point; on the other hand, the latter 
variation produced the best result for (xf). In fact, it is easy to see, by direct integration 
of the second equation ( |4.1| ), that = D/u, which is exactly 0.5 for the values of the 
parameters we chose. By beginning with this value at first order, the methods based on 
perturbation theory had a distinct advantage over the L^L approach — although the values 
of the variational parameters d\ and c?2 were very similar in the MCW and L^L cases. One 
reason for the relatively poor result for (x^) when not using the perturbative schemes, may 
be that the variation in this case is not carried out directly on the quantity of interest. 
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B. The Colored-Noise Example 



Now we move on to consider some of the variational approaches for a more difficult 



situation — the colored-noise problem [15-17]. Consider the following Langevin equation: 



x 



-V'(x) 



(4.6) 



which describes an overdamped particle subject to the force f(x) = —V'(x) and an external 
noise External noise is not intrinsically related to the system's evolution and is typically 
"colored" as opposed to "white," i.e. it is not delta-function correlated. As a consequence, 
many of the techniques and results familiar from the study of Markov processes are not 
applicable. For present purposes, we take £(t) to be Gaussian-distributed with zero mean 
and exponentially correlated: 



mm) = -exp{-\t-t'\M. 

T 



(4.7) 



With this choice, the above one-dimensional non-Markov process ( |4.6|) can be shown to be 
equivalent to the following two-dimensional Markov process |16]: 



± = -V'(x) + £{t), 
1 



T 



T 



(4. 



where rj(t) is Gaussian-distributed with zero mean and this time delta-function correlated: 



{V(t)v(t')) = 2D5(t-t'). 



(4.9) 



The corresponding Fokker-Planck equation is: 
dP(x,£;t) 



dt 



d_ 

dx 



[-V'(x)+Z)P 



+ 



ld_ 



HP 



DdP 



which upon the transformation to the velocity variable £ = x + V'(x) becomes: 



dP(x, X] t) 
dt 



t 2 dx" 2 



V"(x) + 



1 + x 



d_ 

dx 



d 

X TT + 
dx 



V'{x) 



d_ 

dx 



P. 



(4.10) 



(4.11) 
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Again, it is readily confirmed that the potential conditions are not satisfied in this case. 
Nevertheless, there is a solvable case. When the operator L is quadratic (i.e. when the 
drift vector Ai(x) is linear in and the diffusion matrix Bij(x) is independent of the Xj's), one 
can find a closed expression for the stationary distribution even if the potential conditions 
are not met. For the colored- noise problem, this situation occurs if V(x) = uo 2 x 2 /2. The 
associated Fokker-Planck operator is: 



^quad 



D_cP_ 

r 2 dx 2 



+ 



co 2 + 



1 



1+x 



d_ 

dx 



d 

x a~ + 
ox 



v u) 2 x~\ d 



dx' 



(4.12) 



and steady-state distribution is: 
Po oc exp < — 



lj 2 (1 +lu 2 t)^ 2 t(1 + uj 2 t). 2 



2D 



-X 



2D 



x 



(4.13) 



This result can be useful for considering the behavior near the minima for a more general 
potential V(x). 

Hereafter, we take V(x) to be a bistable potential given by V(x) = —x 2 /2 + x 4 /A, with 
the following operator: 



bistable 



D d 2 L n , 1 
+ 3x 2 - 1 + - 
r. 



t 2 dx 2 



. . d 

ox 



X 



d_ 

dx 



x — x 



d_ 

dx 



(4.14) 



For the region immediately surrounding the minima of the bistable V(x) (x = ±1), the 
potential can be described by Taylor expansion truncated at the quadratic order, suggesting 
that 



lim Pq 



exp 



:i + 2r) 
D 



(x - l) 2 



t(1 + 2t) .. 



2D 



-x 



(4.15) 



and similarly around x — — 1. 

As we know of no convenient solvable problem with a bistable potential, let us take the 
L^L approach. All we need is some suitable trial steady-state distribution. We choose a 
variational stationary distribution of the form: 



P(x,x) ~ h(x) exp{ - f(x) - g(x)x 2 }, 



(4.16) 
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so that the integrations over x are readily performed analytically. Note that stability 
requires that g(x) > for all x. Furthermore, the symmetry suggests that f(x), g(x) and 
h(x) are even in x. Calculating the expectation of L^L with respect to P yields: 



{P\L*L\P) 
' (P\P) " 



( I5(g'(x)) 2 3 \ 2Dg{x) 

{64(g(x)) 7/2 Hg(x)) 1/2 [ r 2 
3g'(x) 



- 3x 2 + 1 - 



8G?(*)) 5/2 

1 



4(^)) 3/2 



^(^ +/ / (a:) 2(z 3 -*)<?(*) 



(/i(a;)) 2 exp{-2/(:r)} 



x 



/cfa-^^exp{-2/(x)}] 
J (g(x)) 



(4.17) 



where the i integration has been done. 

We use the following simple polynomials for f(x), g(x) and h(x): 



/(*) = 

g(x) = 
h{x) = 



(l + 2r) 
D 

r(l + 2r) 



[x 



2D 



1 + br 



l) 2 (x 2 
— + ar— 



If 



[x 2 - I) 



CT 



{X 2 - 1) 



(4.18) 



where a, b and c are variational parameters. This motivation for these choices is the follow- 
ing: we know what f(x), g(x) and h(x) should be near x = ±1, and this variation represents 
in some sense the next terms in the Taylor expansions of these functions, expanding simul- 
taneously around x — ±1. (It also recovers the appropriate r — > limit.) If D — 0.1 
and t = 0.6, we find (P\tfL\P) = 0.32015 when a = 0, b = and c = 0, whereas varying 
a, 6 and c yielded a minimum of (P\L^L\P) = 0.037117, an order of magnitude smaller, at 
a = 1.14, b = 2.39 and c = 1.85. The corresponding marginal stationary distribution 

P mar (x) = JdxP(x,x) ~ ^1/2 exp {-f(x) - In /2 + In , (4.19) 

is shown in Figure 1. Note that while there is certainly room for improvement, it does 
capture some of the features such as the shift in the maxima away from x = ±1, due to 
existence of x-dependent prefactors. 
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These results can be compared to those obtained by a systematic small-/} expansion 
using path-integral techniques. The path-integral approach is clearly superior for small D, 
but we expect it to become less reliable as D increases. In contrast, the variational approach 
can, in principle, be applied for any D. 

V. CONCLUSIONS 

In this paper we have discussed a number of variational procedures which may be ap- 
plied to the Fokker-Planck equation. They broadly fall into two classes: those in which 
a Hermitian operator was constructed from the generically non-Hermitian Fokker-Planck 
operator and those which were constructed with reference to an approximation to the full 
problem, which was exactly soluble. Both approaches have advantages and disadvantages. 
For example, in the former approach, where we used the Hermitian operator L^L, twice as 
many operations have to be carried out — the result of first operating with L and then with 
L*. On the other hand, this may be the only method available; in the colored noise example, 
for instance, we were unable to construct a Fokker-Planck operator which was both solvable 
and which contained enough of the essential physics to make it a useful approximation. The 
variants in the second approach were tested on a simple quantum-mechanical example as 
well on a toy Fokker-Planck problem. One of these schemes (MCW) shares with L^L the 
advantage of having a lower bound built in. From the few examples we have investigated, it 
is clear that variational procedures are capable of giving good quantitative, as well as qual- 
itative, results for Fokker-Planck equations. However more work is required to extend the 
range of the problems studied, as well as the order to which the calculations are carried out, 
before definitive statements about the relative usefulness of the various possible approaches 
can be made. 
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APPENDIX: 



In this appendix, some of the calculational details for the toy model are provided. The 
perturbative approaches split the operator L toJ , eq. ( |4.3| ), as follows L toy = L^ y + Ltoy> 

where L[^ y , eq. (O), is the Fokker-Planck operator for the Ornstein-Uhlenbeck process. 

(i) 

toy 



The associated L:ll is 



r(l) 
^toy 



q2 / 
i=l,2 UJj i \ 



_d_ 



(Al) 



It is convenient to work with the transformed operator ^Pq°\x)\ ^ L ^Pq°\x)^J ^ ; ex- 
pressed in terms of annihilation and creation operators, it is: 



-toy = -v(a\ai + a) 2 a 2 ) 



(A2) 



-toy 



( — — u] a\a\ — gd 2 (a 2 + a 2 ) 2 (a\ai + a\a[) 



where X{ = d] (ai + a\) and d/dxi = (a* — a\)/2dy 2 . 
From here one can readily calculate the following: 



(4°VtVtoyl4 0) > 



D 



v - gd 2 



+ 2 



and similarly 



i 

2^2 



D 
d\ 



v - gd 2 



+ 



2u 2 



D 

do 



D 



-i 2 



— V 



+ Wd\ 



V 



1 jNl 

Au 2 ' 



Minimization of (^Vtoy^toyl'^o '') yields 



(0)\ 



(A3) 



(A4) 



(A5) 



while minimization of (^o^l^o^) results in: 



d 1 



2D 2 



and d\ H — d 2 — = 0. 



u + gd 2 ' g< g 

To include the FAC and PMS approaches, we calculate (xf): 



dmdx 2 x, 2 (P (0) (2l) + eP (1) (2l) + e 2 P (2) (2l) + •••)> 



for % — 1, 2. To order e 2 we find: 

(a; 2 ) = di + e 
+e 
= d 2 + e 



v \di 

d i (2g 2 dl , #d 2 £> 



d 2 / -D 



+ 2#d 2 - 



z/di v 



v \d 



— v 



+ [0] • 



(AT) 



(A8) 



(A9) 



FAC gives (x 2 ) = Dvjiy 2 + gD) and (x 2 ,) = D/u at 0(e), and has no physical solution at 
0(e 2 ). PMS, on the other hand, has no physical solution at 0(e) and variation of (x 2 ) at 
0(e 2 ) with respect to d\ and d 2 yields: 

2g 2 d 2 2 gd 2 gD 



+ ^-^ = 



^g 2 did 2 gdi gD _ 

9 9 
V V V 



(A10) 
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Table 1 



Order 


MWC 


FAC 


PMS 


RRL 


1 


0.42957 


0.45428 


0.42927 


0.42927 


2 


0.42242 




(0.42143) 


0.42236 


3 




0.48854 


0.42098 





The results of four variational approaches applied to Eq of the quartic oscillator. 
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Table 2 



Method 


<*?> 


(xl) 


Simulation 


0.265228 


0.500310 


LtL 


0.271875 


0.279693 


MCW 


0.241854 


0.5 


FAC 


0.2 


0.5 


PMS 


0.174307 


0.5 



Variational results for the toy model with D = 0.5, v = 1.0 and g = 3.0. 



23 



FIGURES 



g. 1 The value of ln[P s t(x)] plotted for D = 0.1 and r = 0.6. The dots are from a 
numerical simulation of the Langevin equation and the solid curve is the result of the 
variational calculation. 
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